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Abstract 

This paper discusses the wavelet-based Finite Difference Time Domain (FDTD) method and a tunable high resolution estimator 
with a specific problem of sound wave propagation through phononic crystals. If the band structures of a phononic crystal are 
calculated by the traditional FDTD method combined with the fast Fourier transform (FFT), some disadvantages, such as time 
consuming and the numerical instability of FDTD iterations are encountered. Moreover, good frequency estimation can only be 
ensured by the postprocessing of sufficiently long time series. In this paper, a wavelet-based FDTD and a tunable high 
resolution estimator based on a bank of filters are proposed to overcome these difficulties. Numerical results for two- 
dimensional phononic crystal show that the wavelet-based FDTD method improves the efficiency of the time stepping 
algorithm and the stability of iterations, and tunable high resolution estimator shows the advantages over the FFT-based 
spectral estimation. 
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Introduction 

In 1993, Kushwaha [1] introduced the concept of "phononic crystals" (PNCs), which represented a kind of artificial 
periodic elastic/acoustic structures that exhibited the so-called "phononic band gaps", and the frequency intervals 
where the propagation of acoustic waves (AWs) or elastic waves (EWs) is fully forbidden in all directions. The 
phononic band gaps could be engineered to provide a vibrationless environment for high precision mechanical 
systems in given frequency ranges. Understanding the full band structures is expected to lead to the design of new 
generations of sound shields, filters, transducers, refractive devices such as acoustic lenses and acoustic 
interferometers, etc. Therefore, the theoretical study of PNCs is focused on the band structure calculation. The 
numerical methods developed for phononic energy band calculation include the plane wave expansion (PWE) 
method [2], multiple scattering theory (MST) [3], the wavelet method [4], the finite difference time domain (FDTD) 
method [5] and the finite element method (FEM) [6], etc. Among them, the FDTD method has some noticeable 
advantages, such as its ability to deal with large mismatch systems and its suitability to tackle arbitrarily shaped 
inclusions. 

Despite these advantages, the traditional FDTD method is often time-consuming, poor in convergence and requires 
vast computational resources when simulating phononic crystals with very large acoustic mismatches. Moreover, 
to obtain the frequency bands of a PNC by using the FDTD method, a Fourier transform method must be applied 
into the oversampled FDTD time series. Thus, in order to improve the accuracy of the Fourier transform, a large 
number of FDTD iterations are often needed to generate sufficiently long time series. However, when the 
computation load of each FDTD iteration is relatively big, much computation time is needed to run so many FDTD 
iterations. 

In the past years, several methods combining FDTD algorithms and wavelets have been developed in order to 
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reduce simulation times. Besides, recently, high resolution spectral estimation methods have attracted increasing 
attentions in modern signal processing. These methods can reduce the number of iterations by estimating the 
frequency response of the devices accurately with relatively short FDTD time series. 

In this paper, we present a novel wavelet-based FDTD method combined with a tunable high resolution estimator. 
This scheme combines the favorable characteristics of a FDTD scheme and the possibility to resolve the acoustic 
field on locally different resolutions. Comparisons should be made to clarify the advantages of the new method. 
The aim is to improve the vast computational load and to reduce the number of FDTD iterations by estimating the 
eigenfrequencies of PNCs accurately with relatively short time series, under the circumstance of systems with high 
acoustic impedance ratios. 


Traditional and Wavelet-Based FDTD Method 


An easy way to comply with the journal paper formatting requirements is to use this document as a template and 
simply type your text into it. 


The acoustic wave equations in a 2D phononic system are 
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where r = (x,y); g is compressibility; p is density; p is acoustic pressure; and V - (v x , V y ) is the velocity. Based on 
the Bloch theorem: p = e ,kr U{r,t), v x = e jk T(r ,t) , and v y = e ikr W(r,t) / where k = (k x ,k y ) is the wave vector, and U, V 
and W are periodic functions of the spatial position with the same period as the crystal lattice. Then equations (1) 
and (2) can be rewritten as 
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Traditional FDTD combined with FFT 


In the traditional FFT-based FDTD calculation, the discretized forms of the equations (3)~(5) are as follows 
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where the superscripts (1, m) denote the position (lAx, mAy); l±l/2 and m±l/2 denote position (l±l/2)Ax and 
(m±l/2)Ay respectively; superscripts n and n±l/2 denote time nAt and (n±l/2)At respectively; K' x = {)k x Ax + 2) /(2Ax) 
and K* = (j/cA v f ± 2) /(2Av) ■ When FDTD is used in periodic media, periodic boundary condition must be applied. 

Take the square lattice with period as an example, under this situation the periodic boundary condition requires 
that values of the same field variable should be equal at x=0 and x=a, or y=0 and y=a. After certain initial conditions 
are given, FDTD proceeds with iterations. According to equation (7) (or (8)), the spatial distribution of V (or W) at 
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time (n+l/2)At can be obtained based on the spatial distribution of V (or W) at time (n-l/2)At, and U at time nAt. 
While according to equation (6), the spatial distribution of U at time (n+l)At can be obtained based on the spatial 
distribution of U at time nAt, and the spatial distribution of V (or W) at time (n+l/2)At. The iterations are 
performed step by step. Then the data series Ul,m;n (n=0, 1, 2, . . .) can be obtained for each discrete spatial position 
(lAx, mAy). A disturbance at a random selected position (10Ax, mOAy) can be chosen as the initial condition at the 
start time step (n=0), i.e. Ul,m;0=5(l-10)-5(m-m0). 

V 

To obtain the eigenfrequencies of the phononic crystal corresponding to a given wave vector k , the FDTD time 
series , i.e., the discrete-time series of the value of one field variable on an observational point computed by FDTD 
at each iteration, is usually postprocessed by using the Fourier transform to get its power spectrum. Then the peak 
positions in the power spectrum are identified as the eigenfrequencies. 


Wavelet-based FDTD 

The FDTD can also be derived from the Galerkin method applied in the AW equation with wavelet basis functions. 
We use Biorthogonal wavelet for spatial dependence of the function and Haar wavelet for the time dependence of 
the function, we get algorithms with better properties. The Haar basis function has poor regularity. Therefore, it 
gives rise to oscillations. Using functions that are continuous and have continuous derivatives give better results. 
The Garlerkin method gives the following equation 
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where (I>(z) and <\kz) are the wavelet basis and testing functions. The update equations are as follows 
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Here a modified form of the Courant stability condition should be satisfied 
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The fields are expanded in terms of scaling functions as follows 
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The testing function or dual function used here is the dirac- delta function because of the special structure of PNCs. 
Tunable High Resolution Estimator 

The frequency resolution of the Fourier transform is limited by the finite length of the data series. However, 
improvement of the frequency resolution of the spectral estimation with relatively short time series is an intensely 
discussed topic in modem signal processing. Here, we introduce a parametric high resolution estimator that is 
based on a bank of filters in modern signal processing. The aim of the bank of filters is to process the data in order 
to obtain estimates of the energy spectrum at desired points. Typically, the resulting spectra show significantly 
higher resolution and more robust as compared with traditional linear predictive filtering. The method is usually 
based on the auto-regressive/moving average (ARMA) filters of discrete-time signal {x(n)}. We consider the basic 
problem of estimating its power pectral density as follows 
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where 
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The basic idea of the tunable high-resolution estimator method is to model the given signal as an auto- 
regressive/moving average filters of complexity. Then the power spectrum density of the signal can be calculated 
according to (15). 


Numerical Results and Discussion 

To verify the proposed wavelet-based FDTD method and the frequency postprocessing method, we perform the 
detailed calculations of different PNCs with high wave impedance ratios. Two groups will be compared 
quantitatively: ©The wavelet-based and traditional FDTD; © Tunable high-resolution estimator and the FFT-based 
postprocessing. An intrinsic problem of the FDTD calculation of phononic energy band is that some of the 
eigenfrequencies may be missed due to the limited frequency resolution if FDTD is run with relatively few 
iterations. Accordingly, in this paper the frequency resolution is quantitatively evaluated by the number of missed 
eigenfrequencies in the part of the obtained energy band diagram below a certain normalized frequency. The 
computer programs run on a Dell R900 computing sever (with dual 2.13GFFz CPUs, 64G memory and Linux 
operating system). 

A 2D phononic crystal, which is consisted by the steel square-section columns arranged in a square lattice in the 
host air, is considered in this example. The material parameters are A, ee/ = 1610kglm\ c lsleel =6010 ml s r c tsted = 3150w/ s / 
Pair = ^S / m \ c i,atr =340 ml s . The filling fraction is f=0.5. The band structures for 25*10 grid points are shown in figure 
1. To ensure good convergence of the calculation results, the spatial interval is chosen as Ax=a/150 in this example. 
Fig.l presents the energy band diagrams. The results are obtained by applying the traditional FDTD method (black 
lines in Fig.l) and the present wavelet-based method (black dots in Fig.l). It can be seen that the wavelet-based 
FDTD method gives better results. 



FIG.l THE BAND STRUCTURE OF THE 2D STEEL/AIR PNC OBTAINED BY APPLYING THE TRADITIONAL FDTD METHOD (BLACK 
LINES IN FIG.l) AND THE PRESENT WAVELET-BASED FDTD METHOD (BLACK DOTS IN FIG.l) FOR 25 x 10 GRID POINTS. THE INSET 

IN FIG.l SHOWS THE FIRST BRILLOUIN ZONE OF THE SQUARE LATTICE. 

| 1 flo+T 2 

Moreover, the RMS errors e “ s - (/ j J , 0 S re/ W) dt anc [ iteration time taken are shown in table 1. 

TABLE. 1 RMS ERRORS AND TIME FOR THE TRADITIONAL FDTD AND WAVELET-BASED FDTD 


Method 

Grid 

RMS error 

Time(sec) 

FDTD 

100x40 

0.1102 

7.835 

25x10 

0.3471 

0.757 

Wavelet-based FDTD 

100x40 

0.0451 

10.83 

25x10 

0.2353 

1.021 
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The tunable high resolution estimator and the FFT-based postprocessing will be compared quantitatively. In the 
above examples, 10 observing points are randomly selected to record the FDTD time series. Fig.2 presents the 
energy band diagrams obtained by applying the FFT-based method (black dots in Fig.2a), and the present high 
resolution estimator method (black dots in Fig.2b) with the FDTD time series of length 10,000. The red stars in 
Figs.2a and 2b mark the eigenfrequencies missed by the two postprocessing methods, respectively. So when the 
number of the FDTD iterations is 10,000 which is relatively few, the present postprocessing method has 
significantly higher frequency resolution. 



FIG.2 THE ENERGY BANDS OF THE 2D STEEL/AIR PNC OBTAINED BY APPLYING THE FFT-BASED METHOD (BLACK DOTS IN FIG.2A) 
AND THE TUNABLE HIGH RESOLUTION ESTIMATOR METHOD (BLACK DOTS IN FIG.2B) WITH THE TIME SERIES OF LENGTH 10,000. 
THE RED DOTS MARK THE EIGENFREQUENCIES MISSED BY THE CORRESPONDING POSTPROCESSING METHODS. 

Conclusions 

The advantage of using wavelet-based FDTD over the traditional FDTD is clearly discernible. In the 2D AW system 
we get a more accurate solution for the same number of grid points or we can get the same accuracy with fewer 
number of grid points, thus saving a lot of memory and time. Compared with the classic FFT-based postprocessing 
method, the introduced signal processing method-high resolution estimator can give much better estimation of the 
eigenfrequencies of PNCs when FDTD runs with relatively few FDTD iterations. The combined method of wavelet- 
based FDTD and the proposed postprocessing method may be useful to improve the total efficiency of the inverse 
design and optimization of the AW or EW band gap materials. 
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